As a Matter of State: The Role of Thermodynamics in Magnetohydrodynamic Turbulence

Philipp Grete; Brian W. O’Shea; Kris Beckwith

doi:10.3847/1538-4357/ab5aec

1. Introduction

Magnetic fields are ubiquitous in the universe and have been observed on all scales, from stellar and planetary systems to the intracluster medium. Similarly, many astrophysical systems are expected to be governed by or subject to turbulence simply by the large spatial scales involved (Brandenburg & Lazarian 2013). More generally, magnetized turbulence is thought to play a key role in many astrophysical systems and processes, e.g., magnetic field amplification via the turbulent dynamo (Tobias et al. 2013; Federrath 2016), particle acceleration in shock fronts resulting in cosmic rays (Brunetti & Jones 2015), or the formation of jets (Beckwith et al. 2008) and in accretion disks (Balbus & Hawley 1998).

In the absence of detailed 3D spatiotemporal observations and/or experimental data, numerical simulations are often used to support the interpretation of observations, or in the case of turbulence research, have become one of the major drivers of scientific advances. This pertains, for example, to studying energy dissipation and turbulent energy cascades, which serves to illuminate the physical mechanisms of energy redistribution and the local nature of energy transfer within turbulence (e.g., Yang et al. 2016; Grete et al. 2017; Andrés et al. 2018) or to turbulence modeling (Clark et al. 1979; Germano et al. 1991; Chernyshov et al. 2012; Grete et al. 2016), which allows the incorporation of small-scale turbulent effects and feedback in turbulence simulations that usually are not able to capture the full dynamical range.

From an astrophysical point of view, a range of studies have analyzed turbulence dynamics and statistics in a variety of regimes, with a focus on quantities related to observations enabling inference of statistical properties of the plasma below the observational resolution limit. In particular, in the star formation community significant attention is paid to the relation between density fluctuations and sonic Mach number ( ${M}_{{\rm{s}}}$ ). In the isothermal, supersonic case, the density distribution is well-described by a lognormal distribution and the width of the distribution proportional to the sonic Mach number (Padoan et al. 1997; Passot & Vázquez-Semadeni 1998). Moreover, for a given ${M}_{{\rm{s}}}$ and ${M}_{{\rm{a}}}$ , the standard deviation of the density fluctuations contains information on the effective turbulent production mechanism, with respect to (for example) the effects of different ratios of compressive to rotational modes in the forcing (Federrath et al. 2008). Similarly, the departure from an isothermal equation of state (EOS) has been studied for hydrodynamic turbulence and for a polytropic EOS (Federrath & Banerjee 2015) or an adiabatic EOS (Nolan et al. 2015; Mohapatra & Sharma 2019), indicating additional dependencies, on (for example) the adiabatic index γ, the relation between density fluctuations and ${M}_{{\rm{s}}}$ . In the magnetized case the majority of studies focus on the isothermal case , finding an additional dependency on the ratio of thermal to magnetic pressure ${\beta }_{{\rm{p}}}$ (e.g., Kowal et al. 2007; Padoan & Nordlund 2011; Molina et al. 2012). Overall, a highly dynamic picture of turbulence has been found, challenging our ability to examine the general case.

In this paper, we present adiabatic magnetohydrodynamic simulations of stationary turbulence in the subsonic ( ${M}_{{\rm{s}}}\approx 0.2$ to 0.6), super-Alfvénic ( ${M}_{{\rm{a}}}\approx 1.8$ ), and high ${\beta }_{{\rm{p}}}$ (ratio of thermal to magnetic pressure with $10\lesssim {\beta }_{{\rm{p}}}\lesssim 100$ ) regime. This regime approximates the turbulent intracluster medium (ICM; Brunetti & Lazarian 2007; Brüggen & Vazza 2015) even though our magnetohydrodynamics (MHD) model neglects effects stemming from low collisionality (Schekochihin & Cowley 2006; Schekochihin et al. 2008). However, studies including effects from low collisionality, e.g., through the Chew–Goldberger–Low MHD model, found that they generally introduce only small differences compared to the MHD model, e.g., a small increase in density fluctuations (Kowal et al. 2011; Santos-Lima et al. 2017), but leave imprints on Faraday rotation maps (Nakwacki et al. 2016). Here, we specifically focus on the effects of departure from an isothermal equation of state by varying the adiabatic index γ and the cooling function between simulations.

The rest of this paper is organized as follows. In Section 2, we introduce the numerical setup and the simulations conducted. In Section 3, we present results from analyzing the simulations starting with a high-level overview of the energy spectra, to correlations in the stationary regime, to statistics of distribution functions that can potentially be used to break degeneracies between the thermodynamics and sonic Mach number. The results are discussed in Section 4, and we conclude with a brief outlook in Section 5 on how our findings point the way to future measurements that can be used to better diagnose the properties of turbulence in astrophysical plasmas.

2. Numerical Details

In this work we utilize the equations of compressible, ideal MHD equations

$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\right)=0,\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{\partial }_{t}\rho {\boldsymbol{u}}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\otimes {\boldsymbol{u}}-{\boldsymbol{B}}\otimes {\boldsymbol{B}}+{ \mathcal I }{p}_{\mathrm{tot}}\right)=\rho {\boldsymbol{a}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}-{\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right)=0,\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{\partial }_{t}E+{\rm{\nabla }}\cdot \left(\left(E+{p}_{\mathrm{tot}}\right){\boldsymbol{u}}-{\boldsymbol{B}}\left({\boldsymbol{B}}\cdot {\boldsymbol{u}}\right)\right)=\rho {\boldsymbol{a}}\cdot {\boldsymbol{u}}-{ \mathcal L },\end{eqnarray} \tag{ 4 }$

that are closed by an ideal equation of state ${p}_{\mathrm{th}}=\left(\gamma -1\right)\rho e$ , with γ as the ratio of specific heats. The symbols have their usual meaning, i.e., density ρ, velocity ${\boldsymbol{u}}$ , total pressure ${p}_{\mathrm{tot}}={p}_{\mathrm{th}}+{p}_{{\rm{B}}}$ consisting of thermal pressure ${p}_{\mathrm{th}}$ and magnetic pressure ${p}_{{\rm{B}}}\,=1/2{B}^{2}$ , and magnetic field ${\boldsymbol{B}}$ , which includes a factor $1/\sqrt{4\pi }$ . Cooling is included via ${ \mathcal L }$ and $E=1/2\left(\rho {u}^{2}+{B}^{2}\right)+\rho e$ is the total energy density with specific internal energy e. Vector quantities that are not in boldface refer to the L² norm of the vector and ⨂ denotes the outer product. The details of the acceleration field ${\boldsymbol{a}}$ that we use to mechanically drive our simulations are described in Section 2.3.

2.1. Cooling Functions

The cooling curve for optically thin, astrophysically relevant plasmas is not scale-free, and as such its use is undesirable if we wish to achieve an understanding of non-isothermal turbulence in a broader context. To that end, we use two idealized cooling functions in this work: linear cooling with

$\begin{eqnarray}&&{ \mathcal L }={C}_{\mathrm{cool}}\rho e\propto \rho e,\end{eqnarray} \tag{ 5 }$

and cooling that approximates free–free emission with

$\begin{eqnarray}&&{ \mathcal L }={C}_{\mathrm{cool}}{\rho }^{2}{e}^{1/2}\propto {\rho }^{2}{e}^{1/2}.\end{eqnarray} \tag{ 6 }$

In this idealized setup appropriate units are absorbed in ${C}_{\mathrm{cool}}$ . In the majority of the simulations ${C}_{\mathrm{cool}}$ is chosen to approximately balance turbulent dissipation (which, in turn, balances the energy injection from the forcing). In the case of linear cooling

$\begin{eqnarray}&&\left\langle \rho \right\rangle {U}^{3}/L\approx \left\langle { \mathcal L }\right\rangle \approx {C}_{\mathrm{cool}}\left\langle \rho \right\rangle \left\langle e\right\rangle \approx {C}_{\mathrm{cool}}\left\langle {p}_{\mathrm{th}}\right\rangle /\left(\gamma -1\right),\end{eqnarray} \tag{ 7 }$

where the means, $\left\langle \cdot \right\rangle$ , refer to the spatial mean value in the stationary regime, U is the root mean square velocity in the simulation's stationary regime, and L is the characteristic turbulence length scale.

2.2. Implementation

All simulations were conducted with a modified version of the astrophysical MHD code Athena 4.2 (Stone et al. 2008) using the same numerical scheme consisting of second-order reconstruction with slope-limiting in the primitive variables, an HLLD Riemann solver, constrained transport for the magnetic field, and a MUSCL-Hancock integrator (Stone & Gardiner 2009). Moreover, we used first-order flux correction (Lemaster & Stone 2009; Beckwith & Stone 2011) in cells where the second-order scheme described above results in a negative density or pressure—the integration is repeated using first-order reconstruction. Explicit viscosity and resistivity are not included and thus dissipation is of a numerical nature given the shock capturing finite volume scheme, making the simulations implicit large eddy simulations (Grinstein et al. 2007). Strictly speaking, the equations governing the simulations are not the ideal MHD equations but include implicit dissipative terms. The nature of implicit dissipation in the Athena code was examined by Simon et al. (2009) and Salvesen et al. (2014; see also Beckwith et al. 2019); the work of these authors demonstrated the similarity of these terms to explicit viscosity and resistivity, similar to techniques adopted for inviscid hydrodynamics (Sytine et al. 2000).

In order to achieve a stationary regime with constant Mach number in a driven, adiabatic simulation a mechanism to remove the dissipated energy is required. We implemented a flexible cooling mechanism⁵ approximating optically thin cooling. In addition to the cooling function itself, we added several constraints to the integration cycle.

First, the time step is limited so that the internal energy is not changing by more than 10% per cycle.

Second, a cooling floor is employed in the form of a pressure floor. For all simulations cooling is turned off in a given cell during a given time step if the pressure drops to values of less than 10⁻⁴ in code units (which is ≃10⁻⁴ of the mean initial value in the calculations).

Third, we ported the "entropy fix" of Beckwith & Stone (2011) for relativistic MHD to the non-relativistic case. The entropy fix introduces the entropy as a passive scalar to the set of equations solved. In case the first-order flux correction fails, i.e., if density or thermal pressure are still negative in a cell after a first-order update, the entropy is used to recover positive values. However, this fix was only required in the two simulation that were thermally marginally stable ( ${\mathtt{M}}{\mathtt{0}}.{\mathtt{54}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{58}}}_{\mathrm{Cool}:\mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$ ) and thermally unstable ( ${\mathtt{M}}{\mathtt{0}}.{\mathtt{70}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{37}}}_{\mathrm{Cool}:\mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$ ) (see Section 3.3), and even in those cases only tens out of 1024³ cells were affected per simulation for a very small fraction the time steps.

2.3. Simulations

In total, we conduct 30 simulations. All simulations evolve on a uniform, static, cubic grid with 512³ or 1024³ cells and side length ${L}_{\mathrm{box}}=1$ starting with uniform initial conditions (all in code units) ρ = 1, and ${\boldsymbol{u}}={\boldsymbol{0}}$ . The initial uniform pressure and background magnetic field (in the x-direction and defined via the ratio of thermal to magnetic pressure, ${\beta }_{{\rm{p}}}={p}_{\mathrm{th}}/{p}_{{\rm{B}}}$ ) vary between simulations, as listed in Table 1. A regime of stationary turbulence is reached by a stochastic forcing process that evolves in space and time so that no artificial compressive modes are introduced to the simulation (Grete et al. 2018). The forcing is purely solenoidal, i.e., ${\rm{\nabla }}\cdot {\boldsymbol{a}}=0$ , and the spectrum is parabolic with the peak at k = 2 using normalized wavenumbers (see, e.g., Schmidt et al. 2009, for more details). Thus, the characteristic length is $L=0.5{L}_{\mathrm{box}}$ . Given that all simulations reach a turbulent sonic Mach number of (or close to) ${M}_{{\rm{s}}}=u/{c}_{{\rm{s}}}=0.5$ in the stationary regime with speed of sound ${c}_{{\rm{s}}}=\sqrt{\gamma {p}_{\mathrm{th}}/\rho }$ , we use ${c}_{{\rm{s}}}$ as the characteristic velocity so that the dynamical time $T=1$ . Each simulations is evolved for 10 T, and 10 equally spaced snapshots per dynamical time are stored for analysis. The stationary regime is generally reached after approximately 3T, and we exclude an additional 2T, as a few simulations required more time to reach equilibrium. All statistical results presented below are only covering the stationary regime, i.e., the statistics are calculated over 51 snapshots between $5T\leqslant t\leqslant 10T$ .

Table 1. Overview of the Simulation Parameters

ID	Simulation/Initial Parameters							Stationary Regime
	Resolution	γ	Cooling	${C}_{\mathrm{cool}}$	a	${p}_{\mathrm{th},\mathrm{init}}$	${\beta }_{{\rm{p}},\mathrm{init}}$	${M}_{{\rm{s}}}$	${M}_{{\rm{a}}}$	${\beta }_{{\rm{p}}}$
${\mathtt{M}}{\mathtt{0}}.{\mathtt{23}}{\mathtt{P}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{isoA}}{\mathtt{0}}.{\mathtt{25}}$	512³	1	no	⋯	0.25	1.00	290	0.23(1)	1.93(9)	141(6)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{37}}{\mathtt{P}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{isoA}}{\mathtt{0}}.{\mathtt{56}}$	512³	1	no	⋯	0.56	1.00	71	0.37(1)	1.64(6)	44.7(1.5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{50}}{\mathtt{P}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{isoA}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	1	no	⋯	1.00	1.00	71	0.50(1)	1.59(5)	22.4(1.5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{59}}{\mathtt{P}}{\mathtt{0}}.{\mathtt{74}}{\mathtt{isoA}}{\mathtt{1}}.{\mathtt{00}}$	512³	1	no	⋯	1.00	0.74	53	0.59(2)	1.75(7)	18.7(1.1)

${\mathtt{M}}{\mathtt{0}}.{\mathtt{23}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{73}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{0}}.{\mathtt{26}}$	512³	7/5	free–free	0.025	0.26	0.75	217	0.228(3)	1.77(6)	85.8(2.8)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{23}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{73}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{0}}.{\mathtt{26}}$	512³	7/5	linear	0.018	0.26	0.75	217	0.226(3)	1.77(6)	86(3)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{35}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{72}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{0}}.{\mathtt{51}}$	512³	7/5	free–free	0.067	0.51	0.75	53	0.35(1)	1.66(5)	36.3(2.4)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{35}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{72}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{0}}.{\mathtt{51}}$	512³	7/5	linear	0.049	0.51	0.75	53	0.35(1)	1.67(6)	37.1(2.2)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{50}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{74}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	7/5	free–free	0.165	1.00	0.71	71	0.50(2)	1.80(13)	20.1(5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{50}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{75}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =7/5}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	7/5	linear	0.125	1.00	0.71	71	0.50(2)	1.74(8)	20(1)

${\mathtt{M}}{\mathtt{0}}.{\mathtt{24}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{90}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{39}}$	512³	5/3	free–free	0.051	0.39	0.94	272	0.24(1)	1.98(8)	81(4)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{23}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{72}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{31}}$	512³	5/3	linear	0.039	0.31	0.75	217	0.235(2)	1.92(5)	76.8(1.9)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{35}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{91}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{73}}$	512³	5/3	free–free	0.132	0.73	0.94	67	0.35(1)	1.80(10)	37(3)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{36}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{72}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{61}}$	512³	5/3	linear	0.107	0.61	0.75	53	0.36(1)	1.65(3)	30.3(1.4)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{54}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{58}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	free–free	0.250	1.00	0.60	72	0.54(2)	1.83(6)	16.8(7)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{54}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{57}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	linear	0.300	1.00	0.60	72	0.54(2)	1.73(11)	15.0(1.1)

${\mathtt{M}}{\mathtt{0}}.{\mathtt{37}}{\mathtt{P}}{\mathtt{1}}.{{\mathtt{14}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	free–free	0.200	1.00	1.40	100	0.37(1)	1.89(13)	35.5(1.9)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{36}}{\mathtt{P}}{\mathtt{1}}.{{\mathtt{15}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	free–free	0.200	1.00	1.40	100	0.36(1)	1.72(8)	31.2(1.5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{41}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{91}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	free–free	0.225	1.00	1.20	86	0.41(1)	1.85(7)	27.5(1.3)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{40}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{93}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	free–free	0.225	1.00	1.20	86	0.40(1)	1.75(9)	24.9(1.5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{46}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{73}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	free–free	0.250	1.00	1.00	71	0.46(1)	1.77(10)	20.6(1.6)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{45}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{74}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	free–free	0.250	1.00	1.00	71	0.45(1)	1.7(10)	19.5(1.4)

${\mathtt{M}}{\mathtt{0}}.{\mathtt{70}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{37}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$	1024³	5/3	free–free	0.330	1.00	0.60	43	0.70	1.66	9.81

${\mathtt{M}}{\mathtt{0}}.{\mathtt{53}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{45}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{80}}$	512³	5/3	free–free	0.211	0.80	0.48	34	0.53(2)	1.68(9)	15.0(7)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{48}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{58}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{0}}.{\mathtt{86}}$	512³	5/3	free–free	0.211	0.86	0.60	43	0.48(1)	1.70(4)	18.3(5)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{40}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{93}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	free–free	0.211	1.00	0.94	67	0.40(1)	1.78(3)	26.4(7)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{52}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{59}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	free–free	0.264	1.00	0.60	43	0.52(1)	1.67(3)	15.70(23)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{43}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{97}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{16}}$	512³	5/3	free–free	0.264	1.16	0.94	67	0.43(1)	1.81(7)	23.6(1.1)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{50}}{\mathtt{P}}{\mathtt{1}}.{{\mathtt{01}}}_{\mathrm{Cool}:\ \mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{56}}$	512³	5/3	free–free	0.412	1.56	0.94	67	0.50(1)	1.86(10)	19.7(7)
${\mathtt{M}}{\mathtt{0}}.{\mathtt{46}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{74}}}_{\mathrm{Cool}:\ \mathrm{lin}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}$	512³	5/3	linear	0.220	1.00	0.80	57	0.46(1)	1.74(6)	20.5(7)

Note. The simulations differ by the numerical resolution N³, the adiabatic index γ, the cooling function (no, linear, or free–free) and cooling coefficient C_cool, the root mean square (rms) power in the acceleration field a, the initial thermal pressure ${p}_{\mathrm{th},\mathrm{init}}$ , and the initial ratio of thermal to magnetic pressure ${\beta }_{{\rm{p}},\mathrm{init}}$ . The saturated values of the rms sonic ${M}_{{\rm{s}}}$ and Alfvénic ${M}_{{\rm{a}}}$ Mach number and ${\beta }_{{\rm{p}}}$ are calculated as the mean in the stationary regime between $5T\leqslant t\leqslant 10T$ (with temporal standard deviations in parentheses). The values given for ${\mathtt{M}}{\mathtt{0}}.{\mathtt{70}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{37}}}_{\mathrm{Cool}:\mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$ are the final values at $t=4T$ when runaway cooling is triggered; see Section 3.4 for more details. The simulation IDs M## P## EOS A## are constructed using the sonic Mach number and thermal pressure in the stationary regime, the EOS used, and the forcing amplitude. An ID ending with an H indicates a high-resolution (1024³) simulation versus 512³ without suffix.

Download table as: ASCII Typeset image

In general, the simulations can be separated along different parameter dimensions in order to disentangle different competing effects with respect to parameters. The main parameter dimensions target effects of different thermodynamics (i.e., the EOS and cooling function) and with respect to varying sonic Mach number ${M}_{{\rm{s}}}$ .

With respect to thermodynamic effects, approximately isothermal runs with $\gamma =1.0001$ and no cooling are included as references. These are compared to simulations that employ different cooling functions and ratios of specific heats, γ. More specifically, γ = 5/3 (for a monoatomic gas) and $\gamma =7/5$ (for a diatomic gas) are used, and cooling varies between a linear cooling function and one that approximates free–free emission. This results in sets of five simulations that differ in their thermodynamic properties.

In addition, these sets are compared at different (sub)sonic Mach numbers, which is realized by either varying u through varying forcing amplitudes and/or by varying the mean thermal pressure. An overview of all simulations and their parameters is given in Table 1.

3. Results

3.1. Overview

All simulations reach an approximately stationary state in the subsonic ( $0.2\lesssim {M}_{{\rm{s}}}\lesssim 0.6$ ), super-Alfvénic ( ${M}_{{\rm{a}}}=u/{v}_{{\rm{A}}}\approx 1.8$ with Alfvén velocity ${v}_{{\rm{A}}}=B/\sqrt{\rho }$ ), high ${\beta }_{{\rm{p}}}$ ( $10\lesssim {\beta }_{{\rm{p}}}\lesssim 100$ ) regime. In other words, using these dimensionless numbers as proxies means that the kinetic energy, on average, is slightly lower than the thermal energy and slightly larger than the magnetic energy, and that the thermal energy (or pressure) is larger than the magnetic energy (or pressure).

For reference, the temporal evolution of ${M}_{{\rm{s}}}$ , ${M}_{{\rm{a}}}$ , and ${\beta }_{{\rm{p}}}$ for five simulations with ${M}_{{\rm{s}}}\approx 0.5$ but with varying cooling function and adiabatic index γ is illustrated in Figure 1. The transient phase in which the initial conditions evolve toward the stationary regime under constant driving lasts for about three dynamical times (3 T). In general, all data in the stationary regime presented in the following spans the temporal mean (and variations) between $5T\leqslant t\leqslant 10T$ , which excludes an additional $2T$ between $3T\leqslant t\leqslant 5T$ , as few simulations took longer to reach approximate equilibria.

**Figure 1.** Temporal evolution of the sonic Mach number ${M}_{{\rm{s}}}$ , the Alfvènic Mach number ${M}_{{\rm{a}}}$ , and the ratio of thermal to magnetic pressure ${\beta }_{{\rm{p}}}$ for five simulations with ${M}_{{\rm{s}}}\approx 0.5$ in the stationary regime. The simulations differ by the cooling function (none, linear, and free–free) used and the adiabatic index γ.
Download figure:
Standard image High-resolution image

Despite varying thermodynamics (isothermal EOS, adiabatic EOS with γ = 5/3 and $\gamma =7/5$ , and linear and free–free cooling) the five simulations in Figure 1 reach practically identical ${M}_{{\rm{s}}}$ , ${M}_{{\rm{a}}}$ , and ${\beta }_{{\rm{p}}}$ in the stationary regime (see Table 1).

Similarly, the mean kinetic, magnetic, and internal energy spectra⁶ are also identical as shown in Figure 2 (left column) for the same five simulations. The kinetic energy spectrum exhibits a power-law scaling within the wavenumber range 4 ≲ k ≲ 40. No clear power-law scaling is observed in the magnetic and internal energy spectra. The right column of Figure 2 shows the energy spectra for five simulations with free–free cooling and γ = 5/3 but with varying $0.24\lesssim {M}_{{\rm{s}}}\lesssim 0.54$ . Again, all spectra are identical between the simulations apart from the shorter extent at high wavenumbers of the simulation run at 512³ compared to the others run at 1024³. General differences in the raw power (i.e., vertical offsets due to different numerical values of, e.g., u in the simulations) have been removed by normalizing the area under the spectra to unity. This emphasizes the identical shape of the power spectra of all simulations.

**Figure 2.** Mean energy spectra of kinetic energy (top row), magnetic energy (middle row), and internal energy (bottom row). The mean is taken over the stationary regime $t\gt 5T$ and the spectra are compensated by power laws of exponents 4/3, 5/3, and 4/3, respectively. The kinetic energy spectrum is calculated based on the Fourier transform of $\sqrt{\rho }u$ and the internal energy spectrum based on $\sqrt{\rho }{c}_{s}$ . All spectra are normalized to unit area under the curve. The left column shows the same simulations as in Figure 1, i.e., ${M}_{{\rm{s}}}\approx 0.5$ with different cooling functions, and the right column shows only simulations with free–free cooling and γ = 5/3 but with varying ${M}_{{\rm{s}}}$ . Simulations with labels ending in `H` were run at 1024³ cell resolution, with all other calculations at 512³.
Download figure:
Standard image High-resolution image

3.2. Correlations

Similar to temporal evolution of the Mach numbers, the correlation coefficient between thermal ( ${p}_{\mathrm{th}}$ ) and magnetic pressure ( ${p}_{{\rm{B}}}$ ) for the five simulations with ${M}_{{\rm{s}}}\approx 0.5$ and varying EOS and cooling settles to the same value of ≈−0.8 in the stationary regime (see the bottom left panel in Figure 3). In contrast to this, different EOS and cooling functions result in different correlation coefficients between the density field (ρ) and the magnetic field strength (B), as illustrated in the top left panel of Figure 3. The isothermal reference case exhibits a strong anticorrelation of −0.81(2), as previously observed in similar simulations (Yoon et al. 2016; Grete et al. 2018). The anticorrelation is weakened when departing from an isothermal EOS. For γ = 7/5 the coefficient is ≈−0.71 independent of the cooling function, and for γ = 5/3 it is −0.64(3) in the case of linear cooling and −0.55(3) for free–free cooling.

**Figure 3.** Correlation between the density field ρ and the magnetic field strength B (top) and the correlation between thermal and magnetic pressure (bottom). The left column shows the temporal evolution of the correlations for the same five simulations as in Figure 1, i.e., simulations with ${M}_{{\rm{s}}}\approx 0.5$ but with varying EOS and cooling. For the purpose of illustration, the dotted line in the top left panel indicates the mean value and the shaded area the standard deviation of that quantity over time ( $5T\leqslant t\leqslant 10T$ ) as it is used in the right column. The right column shows the mean correlation coefficients vs. sonic Mach number ${M}_{{\rm{s}}}$ in the stationary regime. Each data point corresponds to 1 of 30 simulations total. The horizontal and vertical lines for each symbol (usually within the bounds of the symbol) correspond to standard deviation over time. Filled symbols correspond to simulations run at a resolution of 1024³ and empty symbols to a resolution of 512³. The gray area highlights several simulation pairs running with identical parameters at different resolution that are used for convergence analysis (see the Appendix).
Download figure:
Standard image High-resolution image

These trends are also observed for different sonic Mach numbers, as shown in the right column of Figure 3. Here, the mean and standard deviation of the ρ–B and ${p}_{\mathrm{th}}$ – ${p}_{{\rm{B}}}$ correlation coefficients in the stationary regime are illustrated versus sonic Mach number for all 30 simulations. In the regime presented, the ${p}_{\mathrm{th}}$ – ${p}_{{\rm{B}}}$ correlation coefficient (≈−0.8) is practically independent of sonic Mach number, EOS, and cooling, with a very weak trend toward weaker anticorrelation with increasing Mach number. Overall, the thermal and magnetic pressures are highly anticorrelated. This indicates a total pressure equilibrium (see also ${p}_{\mathrm{tot}}$ distributions in the following Section 3.3).

The individual ρ–B correlation coefficients are predominately determined by the EOS (here, via γ) and the cooling function used; see the top right panel of Figure 3. A higher adiabatic index ( $\gamma =1.0001\to 7/5\to 5/3$ ) results in weaker anticorrelations. Moreover, free–free cooling results in slightly weaker anticorrelations compared to linear cooling, but this effect is mostly visible in the γ = 5/3 simulations. For example, for ${M}_{{\rm{s}}}\approx 0.35$ the ρ–B correlation coefficients are −0.84 in the isothermal case, −0.80 in the γ = 7/5 case with linear cooling and −0.78 with free–free cooling, and −0.73 in the γ = 5/3 case with linear cooling and −0.69 with free–free cooling. Finally, there is an indication that higher numerical resolution also results in slightly weaker ρ–B anticorrelations (of about the same order as the cooling functions). The Appendix presents a discussion of these results with simulation resolution.

3.3. Probability Density Functions (PDFs)

Similar to the ρ–B correlation coefficients, different γ and different cooling functions lead to systematically changing statistics in other quantities. Figure 4 shows the mean PDFs of the density $\mathrm{ln}\rho$ , the normalized thermal pressure ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle$ , the normalized total pressure ${p}_{{\rm{t}}{\rm{o}}{\rm{t}}}/ \langle {p}_{{\rm{t}}{\rm{o}}{\rm{t}}} \rangle$ , and normalized deviation of the derived line-of-sight (LOS) magnetic field strength to the actual one $({B}_{\mathrm{LOS}}-{B}_{0})/{B}_{0}$ . The latter is derived from rotation measures⁷ via

$\begin{eqnarray}&&{B}_{\mathrm{LOS}}=\displaystyle \frac{{\int }_{L}\rho \left(l\right){B}_{\parallel }\left(l\right){\rm{d}}l}{{\int }_{L}\rho \left(l\right){\rm{d}}l},\end{eqnarray} \tag{ 8 }$

with ${B}_{\parallel }$ being the LOS component of the magnetic field.

$\mathrm{ln}\rho $ — **Figure 4.** Mean probability density functions (PDF) of the density $\mathrm{ln}\rho$ , the normalized thermal pressure ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle$ , the normalized total pressure ${p}_{\mathrm{tot}}/\left\langle {p}_{\mathrm{tot}}\right\rangle$ , and normalized deviation of the derived line-of-sight magnetic field strength to the actual one $({B}_{\mathrm{LOS}}-{B}_{0})/{B}_{0}$ . Normalization is applied with respect to the mean value. The mean is taken over the stationary regime ( $t\gt 5T$ ) and the shaded regions indicate the standard deviation of the PDFs over time. The rows show the same simulation as in Figure 1, i.e., ${M}_{{\rm{s}}}\approx 0.5$ with different cooling functions. The bottom row depicts only simulations with free–free cooling and γ = 5/3 but with varying resolution, ${M}_{{\rm{s}}}$ , and power in the forcing field.
Download figure:
Standard image High-resolution image

The top row in Figure 4 shows five simulations with ${M}_{{\rm{s}}}\approx 0.5$ and with different EOS and cooling functions. Differences in the shape of the PDFs of the density, thermal pressure and total pressure between the isothermal case, γ = 7/5, and γ = 5/3 are immediately apparent. For example, with higher γ the PDF of $\mathrm{ln}\rho$ becomes less skewed and all three PDFs become broader. Differences between linear and free–free cooling are more subtle as discussed below. No clear signal between different EOS is observed in the derived LOS magnetic field strengths.

The bottom row in Figure 4 shows five simulation with γ = 5/3 and free–free cooling but with varying sonic Mach number ( $0.25\lesssim {M}_{{\rm{s}}}\lesssim 0.55$ ). Again, differences in the shapes of the PDFs of the density, thermal pressure, and total pressure are apparent. With increasing sonic Mach number all PDFs become broader. In addition, the PDF of the thermal pressure becomes less skewed with increasing ${M}_{{\rm{s}}}$ , while the PDF of the total pressure remains mostly symmetric. In fact, these differences with ${M}_{{\rm{s}}}$ are much more pronounced (see the scaling of the y-axis), suggesting that the sonic Mach number is the dominant effect compared to changes in the EOS and cooling. Again, no significant differences in the PDFs of the LOS magnetic field strength are observed.

In order to further quantify the results, we calculate the statistical moments (mean, standard deviation, skewness, and kurtosis) of these four quantities for all snapshots of all 30 simulations. The skewness is calculated as

$\begin{eqnarray}&&\mathrm{skew}x=\displaystyle \frac{\left\langle {\left(x-\left\langle x\right\rangle \right)}^{3}\right\rangle }{{\sigma }^{3}\left(x\right)}\end{eqnarray} \tag{ 9 }$

with standard deviation σ, and the (Fisher) kurtosis is calculated as

$\begin{eqnarray}&&\mathrm{kurt}x=\displaystyle \frac{\left\langle {\left(x-\left\langle x\right\rangle \right)}^{4}\right\rangle }{{\sigma }^{4}\left(x\right)}-3.\end{eqnarray} \tag{ 10 }$

The mean (over time) statistical moments, including their standard deviations (over time) versus sonic Mach number ${M}_{{\rm{s}}}$ are illustrated in Figure 5. Moreover, in cases where absolute correlation between a quantity x and ${M}_{{\rm{s}}}$ is larger than 0.9, we perform a linear fit with

$\begin{eqnarray}&&x={{mM}}_{{\rm{s}}}+b.\end{eqnarray} \tag{ 11 }$

The regressions are done over all outputs of all simulations employing a particular combination of EOS and cooling, e.g., for γ = 7/5 with linear cooling $3\times 51=153$ data points are taken into account or $4\times 51=204$ in the isothermal case. The slope m of the fit and the correlation coefficient are given in Table 2. Note that given the limited range of ${M}_{{\rm{s}}}$ these fits need to be interpreted with care and we primarily use them here in order to quantify differences (or the absence thereof) between different equations of state and cooling functions.

Table 2. Overview of the Numerical Values of the Linear Fitting Results Illustrated in Figure 5

Quan.	Stat.		Isoth.	γ = 7/5 lin.	γ = 7/5 ff.	$\gamma =5/3$ lin.	$\gamma =5/3$ ff.
$\mathrm{ln}\rho /\left\langle \rho \right\rangle$	mean	m	−0.063(1)	−0.052(1)	−0.060(1)	−0.073(1)	−0.103(1)
		Corr	−0.98	−0.98	−0.97	−0.97	−0.95

$\mathrm{ln}\rho /\left\langle \rho \right\rangle$	std.	m	0.483(7)	0.495(5)	0.536(6)	0.575(6)	0.664(5)
		Corr	0.99	0.99	0.99	0.99	0.98

${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle$	std.	m	0.420(6)	0.523(4)	0.528(5)	0.624(5)	0.607(4)
		Corr	0.99	0.99	0.99	0.99	0.98

${p}_{\mathrm{th}}$	skew	m	2.57(10)	2.07(7)	2.18(7)	2.27(5)	2.68(5)
		Corr	0.90	0.92	0.94	0.95	0.89

${p}_{\mathrm{th}}$	kurt.	m	−7.4(4)	−8.1(3)	−8.4(3)	−8.1(2)	−8.7(2)
		Corr	0.86	0.93	0.94	0.95	0.91

${p}_{\mathrm{tot}}/\left\langle {p}_{\mathrm{tot}}\right\rangle$	std.	m	0.250(4)	0.252(2)	0.262(3)	0.289(4)	0.288(3)
		Corr	0.98	0.99	0.99	0.99	0.97

Note. m is the slope of the linear fit (including standard deviation) and Corr is the correlation coefficient.

Download table as: ASCII Typeset image

Both the mean and the standard deviation ${\sigma }_{\mathrm{ln}\rho }$ exhibit a high (anti)correlation with ${M}_{{\rm{s}}}$ of ≤−0.95 and ≥0.98, respectively. The trend of broader $\mathrm{ln}\rho$ distributions with increasing ${M}_{{\rm{s}}}$ observed in Figure 4 holds across all combinations of EOS and cooling. Based on the slopes, this trend is more pronounced both with larger γ and with cooling that is more sensitive to ρ. In the isothermal reference case the slope is shallower (−0.48) compared to γ = 7/5, with linear (−0.50) and free–free (−0.54) cooling, and to γ = 5/3 with linear (−0.58) and free–free (0.66) cooling. For the skewness and kurtosis no dependency on ${M}_{{\rm{s}}}$ is observed. However, the distributions generally separate for different γ and become less skewed and less broad with increasing γ.

All statistical moments of the normalized thermal pressure (second row in Figure 5) vary with ${M}_{{\rm{s}}}$ . Similar to the density distributions, the standard deviation of the thermal pressure tightly depends on ${M}_{{\rm{s}}}$ (correlation coefficient ≥0.98) and exhibits an additional (weaker) dependency on γ. With increasing γ the slopes are getting steeper, from ≈0.42 in the isothermal case, to ≈0.52 for γ = 7/5, to ≈0.61 for γ = 5/3 (with no pronounced difference between cooling functions). For the skewness and kurtosis the correlations with ${M}_{{\rm{s}}}$ are generally weaker but still pronounced (≥0.86) and a clear trend differentiating EOSs and cooling functions is not observed.

In contrast to this, the skewness and kurtosis of the normalized total pressure distributions (third row in Figure 5) are independent of ${M}_{{\rm{s}}}$ and also independent of γ or a cooling function. However, the standard deviation is again tightly correlated (≥0.97) with ${M}_{{\rm{s}}}$ , and similarly to the thermal pressure and density, shows an additional (weaker) dependency on γ with steeper slopes for larger γ.

Finally, the statistical moments of the normalized derived LOS magnetic field strength distributions (bottom row in Figure 5) generally exhibit no clear trend with ${M}_{{\rm{s}}}$ , EOS, or cooling. A weak trend is seen only in the mean value for a stronger underestimation of the field strength with increasing ${M}_{{\rm{s}}}$ , but the scatter is too large to large to make a definite statement.

3.4. Pressure–Density Dynamics and Thermal Stability

In order to first understand the individual distributions presented in the previous section, the mean 2D PDFs of thermal pressure versus density are illustrated in Figure 6.

**Figure 6.** Mean 2D PDFs of the normalized thermal pressure ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle$ vs. density. The mean is covering 51 snapshots in the stationary regime $5T\leqslant t\leqslant 10T$ . For reference, the gray dashed lines indicates ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle =\rho$ . The top row illustrates the same simulations as in Figure 1, i.e., ${M}_{{\rm{s}}}\approx 0.5$ with different cooling functions, and the bottom row shows simulations with free–free cooling and γ = 5/3 but varying ${M}_{{\rm{s}}}$ (through lowering the mean pressure from left to right). Instead of the mean PDF the bottom right figure shows the final PDF at $t=4T$ when runaway cooling is triggered; see Section 3.4 for more details.
Download figure:
Standard image High-resolution image

${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle $ — **Figure 6.** Mean 2D PDFs of the normalized thermal pressure ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle$ vs. density. The mean is covering 51 snapshots in the stationary regime $5T\leqslant t\leqslant 10T$ . For reference, the gray dashed lines indicates ${p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle =\rho$ . The top row illustrates the same simulations as in Figure 1, i.e., ${M}_{{\rm{s}}}\approx 0.5$ with different cooling functions, and the bottom row shows simulations with free–free cooling and γ = 5/3 but varying ${M}_{{\rm{s}}}$ (through lowering the mean pressure from left to right). Instead of the mean PDF the bottom right figure shows the final PDF at $t=4T$ when runaway cooling is triggered; see Section 3.4 for more details.
Download figure:
Standard image High-resolution image

The top row depicts the PDFs for the simulations at ${M}_{{\rm{s}}}\approx 0.5$ with varying γ and cooling. In general, the distributions are extended around the isothermal reference line ( $p\propto \rho$ ) as expected given the chosen balance between turbulent dissipation and cooling. With higher γ the distributions are getting broader in both dimensions. Moreover, free–free cooling leads to an additional broadening in the density dimension in both cases for γ = 7/5 and $\gamma =5/3$ . The most extreme case (γ = 5/3 with free–free cooling in the top right panel) exhibits broad density tails, as the simulation is thermally marginally stable.

To further illustrate the transition to a thermally unstable regime the bottom row in Figure 6 shows the mean 2D PDFs only for simulations with γ = 5/3 and free–free cooling but with increasing ${M}_{{\rm{s}}}$ (going from 0.36 to $\approx 0.7$ ). The increasing ${M}_{{\rm{s}}}$ (for the same forcing amplitude) is achieved by lowering the mean thermal pressure in the simulations. With increasing sonic Mach number the distributions are getting broader in both dimensions. This can be attributed to the increasing width of the density PDF with ${M}_{{\rm{s}}}$ (see Section 3.3), which is enabled by decreasing pressure support against compression. The bottom right panel shows simulation ${\mathtt{M}}{\mathtt{0}}.{\mathtt{70}}{\mathtt{P}}{\mathtt{0}}.{{\mathtt{37}}}_{\mathrm{Cool}:\mathrm{ff}}^{\gamma =5/3}{\mathtt{A}}{\mathtt{1}}.{\mathtt{00}}{\mathtt{H}}$ for which the pressure support is insufficient to prevent runaway cooling resulting in extended high-density tails.

4. Discussion

4.1. Correlations and Relevance to Observations

The correlation between the density ρ and magnetic field strength B is astrophysically relevant to the LOS magnetic field strength measurement via Faraday rotation. Only for uncorrelated fields and in the isothermal case the derived strength is exact (Beck & Wielebinski 2013). Here, we observe that the ρ–B correlation depends on the adiabatic index γ, i.e., there is a clear departure from the isothermal case. The strong anticorrelation observed in isothermal simulation weakens with larger γ. For isothermal simulations it was additionally observed that the correlation depends on the sonic Mach number ${M}_{{\rm{s}}}$ (especially when going to the supersonic regime) and the correlation time of the forcing (Yoon et al. 2016; Grete et al. 2018; Beckwith et al. 2019). However, the mean deviation of the derived LOS magnetic field from the exact one in our simulations is at most 10%, with a trend of the deviation becoming more significant going from ${M}_{{\rm{s}}}\approx 0.2$ to ≈0.6 independent of different thermodynamics. Thus, the resulting deviation for ICM-like plasmas is likely below the observational uncertainties.

Overall, ρ–B are anticorrelated across all sonic Mach numbers presented. For isothermal MHD Passot & Vázquez-Semadeni (2003) showed that this anticorrelation is indicative of dynamics governed by slow magnetosonic modes. The dependency on γ in the ρ–B correlation observed in this paper suggests that there exists a richer mix of modes when departing from an isothermal EOS.

In contrast to the ρ–B correlation, the correlation between thermal and magnetic pressure is independent of γ and cooling, i.e., the correlation coefficient of the isothermal simulation is indistinguishable from the adiabatic simulations with cooling. This suggests that all simulations are governed by a total pressure equilibrium, ${p}_{\mathrm{th}}+{p}_{{\rm{B}}}={p}_{\mathrm{tot}}\approx \mathrm{const}.$ .

4.2. ${\sigma }_{\rho }$ – ${M}_{{\rm{s}}}$ Relation and Comparison to Previous Work

The presented work covers isothermal and non-isothermal magnetized stationary turbulence with varying γ and cooling functions over a range of ${M}_{{\rm{s}}}$ in the subsonic regime. The majority of previous related work comes from the star formation community and targets isothermal, (magneto)hydrodynamic, supersonic turbulence.

Initial work on the relation between density variations and the sonic Mach number goes back to Padoan et al. (1997) and Passot & Vázquez-Semadeni (1998), who derived and tested numerically the linear relation

$\begin{eqnarray}&&{\sigma }_{\rho /\left\langle \rho \right\rangle }={{bM}}_{{\rm{s}}}.\end{eqnarray} \tag{ 12 }$

In the case of isothermal hydrodynamic turbulence, Federrath et al. (2008) later showed that the proportionality constant b varies depending on the modes employed in the forcing (between ≈0.3 for purely solenoidal forcing and ≈1 for purely compressive forcing).

Qualitatively, we also find a linear relation between density variations⁸ and ${M}_{{\rm{s}}}$ . Moreover, we find that the slope of the linear relation depends on both the adiabatic index γ (steeper with larger γ) and the cooling employed for identical, purely solenoidal forcing. This suggests that departure from an isothermal regime adds additional complexity to the relation, which is also found by Nolan et al. (2015) in the hydrodynamic case. The latter presents both numerical results and a theoretical model that predicts steeper slopes for larger γ in the subsonic regime.

In the MHD case adjustments to the relation have been reported (by, e.g., Padoan & Nordlund 2011; Molina et al. 2012) that take the ratio of thermal to magnetic pressure, ${\beta }_{{\rm{p}}}$ , into account. However, the adjustment is of the order of $1/\sqrt{1+{\beta }_{{\rm{p}}}^{-1}}$ . Given that $10\lesssim {\beta }_{{\rm{p}}}\lesssim 100$ in all of our simulations, this correction would contribute at most a few percent to our results and is thus negligible.

Kowal et al. (2007) studied density fluctuations in isothermal MHD turbulence, including higher-order statistical moments such as skewness and kurtosis. However, the random (uncorrelated) forcing employed in Kowal et al. (2007) leads to the excitation of compressive modes (despite solenoidal forcing) in the subsonic regime that systematically affect several statistics including the correlation between density and magnetic field strength or the density PDF (Yoon et al. 2016; Grete et al. 2018). This systematic effect renders a direct comparison difficult as the results are dependent on multiple parameters. For example, a shorter autocorrelation in the forcing requires a larger forcing amplitude to reach the same sonic Mach number, resulting in more power in compressive modes. In turn, this changes the density PDF for identical sonic Mach numbers, and thus is complementary to the changes described for varying Mach number and EOS in this manuscript.

Finally, it should be noted that our results are not in agreement with recently published results by Mohapatra & Sharma (2019), who conducted adiabatic hydrodynamic simulations with and without heating/cooling. They find ${\sigma }_{{p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle }\propto {\sigma }_{\rho /\left\langle \rho \right\rangle }\propto {M}_{{\rm{s}}}^{2}$ in the subsonic regime without cooling, and ${\sigma }_{{p}_{\mathrm{th}}/\left\langle {p}_{\mathrm{th}}\right\rangle }\propto {M}_{{\rm{s}}}^{2}$ and ${\sigma }_{\rho /\left\langle \rho \right\rangle }\gt {M}_{{\rm{s}}}$ with cooling, whereas we find linear relationships for both density and pressure fluctuations. Given the differences in the setup, e.g., MHD versus HD, idealized cooling versus realistic cooling curve, thermally unstable versus stable, and heating only via turbulent dissipation versus turbulent dissipation and explicit heating, the observed differences in the results may stem from a variety of sources or a combination thereof. As a result, we refrain from a more detailed, purely speculative comparison between the results presented here and those of Mohapatra & Sharma (2019).

4.3. Limitations

Given the idealized nature of this work several items need to be kept in mind when interpreting or extrapolating from the results.

This pertains, for example, to the idealized cooling functions that in their current form only approximate subregimes of a realistic cooling function. Similarly, given the monotonic shape of the cooling functions and the targeted balance between turbulent dissipation and cooling (to achieve stationary turbulence) prevents the development of multi-phase flows. Thus, the results presented provide a qualitative view on the effects of different cooling functions and EOS. For detailed predictions in specific environments such as different phases in the ISM more realistic cooling functions should be employed.

In addition, the sampled parameter space is mostly targeted at ICM-like regimes, i.e., subsonic, super-Alfvénic, and high ${\beta }_{{\rm{p}}}$ turbulence, though neglecting effects from low collisionality in the ICM that can also alter statistical moments of, for example, the density distribution (Schekochihin & Cowley 2006; Kowal et al. 2011). While several clear trends in the ρ–B correlations with varying thermodynamics and in the distribution functions of ρ, ${p}_{\mathrm{th}}$ , and ${p}_{{\rm{B}}}$ with varying ${M}_{{\rm{s}}}$ have been observed, the resulting relations should be handled with care—especially in extrapolating to the supersonic regime.

5. Conclusions

In this paper, we systematically studied how the departure from an isothermal EOS affects stationary magnetohydrodynamic turbulence. In total, we conducted 30 numerical simulations with varying adiabatic index γ with $\gamma \to 1$ for an approximately isothermal gas as reference case, γ = 7/5 for a diatomic gas, and $\gamma =5/3$ for a monoatomic gas. Moreover, we employed two idealized cooling function (linear cooling with $\dot{E}\propto \rho e$ and approximate free–free emission with $\dot{E}\,\propto {\rho }^{2}\sqrt{e}$ ) in order to maintain stationary turbulence with a constant Mach number. All simulations are subsonic ( ${M}_{{\rm{s}}}\approx 0.2$ to 0.6), super-Alfvénic ( ${M}_{{\rm{a}}}\approx 1.8$ ), and high ${\beta }_{{\rm{p}}}$ (ratio of thermal to magnetic pressure with $10\lesssim {\beta }_{{\rm{p}}}\lesssim 100$ )—a regime found, for example, in the intracluster medium.

In this regime, we find that the kinetic, magnetic, and internal energy spectra are practically unaffected by the thermodynamics and the sonic Mach number (apart from the normalization). Moreover, the thermal and magnetic pressures are strongly anticorrelated (correlation coefficient ≲−0.8) independent of γ and cooling, and only exhibit a weak trend toward weaker anticorrelation with increasing ${M}_{{\rm{s}}}$ . In contrast to this, the correlation between density ρ and magnetic field strength, B (which, again, are anticorrelated) shows a dependency on γ. The correlation coefficient of ≈−0.8 in the isothermal reference case gets weaker with larger γ up to ≈−0.55 for γ = 5/3 with free–free cooling. Departing from an isothermal EOS allows independent thermal pressure and density variations. Thus, for a fixed, strong ${p}_{\mathrm{th}}$ – ${p}_{{\rm{B}}}$ anticorrelation (associated with a total pressure equilibrium) an adiabatic EOS naturally reduces the ρ–B correlation by construction.

Similarly, we find dependencies on γ in multiple distribution functions. However, these dependencies are typically subdominant with respect to the overall trend with ${M}_{{\rm{s}}}$ . For example, we find linear relations for an increase of density fluctuations, thermal and total pressure fluctuations, and the skewness of the thermal pressure distribution with increasing ${M}_{{\rm{s}}}$ . A larger γ and a cooling function with stronger density dependency generally result in slightly steeper slopes of these linear relations.

Overall, this results in degeneracy in inferring, for example, Mach numbers from observed distributions without knowing the governing thermodynamics of the observed system. However, we suggest that higher-order statistics (e.g., the skewness of the density distribution) are less dependent on ${M}_{{\rm{s}}}$ and predominantly determined by γ. Thus, there is hope that this degeneracy can be resolved. To do so would require simulations that span a substantially broader and more complex parameter space, which we leave to future work.

The authors thank Wolfram Schmidt and Christoph Federrath for useful discussions. P.G. and B.W.O. acknowledge funding by NASA Astrophysics Theory Program grant #NNX15AP39G. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. B.W.O. acknowledges additional funding by NSF AAG grant #1514700. The simulations were run on the NASA Pleiades supercomputer through allocation SMD-16-7720 and on the Comet supercomputer as part of the Extreme Science and Engineering Discovery Environment (XSEDE Towns et al. 2014), which is supported by National Science Foundation grant number ACI-1548562, through allocation #TG-AST090040. Athena is developed by a large number of independent researchers from numerous institutions around the world. Their commitment to open science has helped make this work possible. SAND Number: SAND2019-14747 J.

Software: Modified version of Athena (v4.2 Stone et al. 2008) available at https://github.com/pgrete/Athena-Cversion (3a7c300). Matplotlib (Hunter 2007). NumPy (van der Walt et al. 2011). mpi4py (Dalcín et al. 2005). mpi4py-fft (Dalcin et al. 2019).

Appendix: Convergence of Simulations

All simulations presented in this paper were conducted at a grid resolution of either 512³ or 1024³ grid cells, with the latter indicated by a ${\mathtt{H}}$ suffix in the simulation ID. While small differences between simulations with identical parameters but different resolutions were observed in the correlations (see Section 3.2), the dominant effect determining the statistics discussed in the paper is related to varying sonic Mach number. Moreover, the PDFs are converged with resolution, as illustrated in Figure 7, where three sets of simulations with identical ${M}_{{\rm{s}}}$ ( $\approx 0.36,0.4,\mathrm{and}\,0.45$ ) for both resolutions are shown. In general, the PDFs for the same ${M}_{{\rm{s}}}$ (i.e., same color) are on top of each other, i.e., the solid lines for simulations at 1024³ are below the dotted lines of simulations at 512³, illustrating convergence.